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Abstract. Equipping the probability space with a local Dirichlet form with square 
field operator T and generator A allows to improve Monte Carlo computations of 
expectations, densities, and conditional expectations, as soon as we are able to 
simulate a random variable X together with T[X] and We give examples on 

the Wiener space, on the Poisson space and on the Monte Carlo space. When X 
is real-valued we give an explicit formula yielding the density at the speed of the 
law of large numbers. 
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1. INTRODUCTION 

Dirichlet forms techniques have shown their efficiency in order to obtain ex- 
istence of densities under weak hypotheses, especially for stochastic differential 
equations with Lipschitz coefficients (cf [6]). We show here their utility for speed- 
ing up Monte Carlo methods especially for the computation of such densities. 

In the whole article the framework is an error structure (O, A, IP, ID, T), i.e. a 
probability space equipped with a local Dirichlet form with square field operator 
(cf [7], [4]). We denote £ the associated Dirichlet form given by £[u] = \ JT[u] dP 
and (A, T>A) the generator linked with £ by the relation £[u,v] = — < A[u],v > 
\Ju G VA,Vv G ID. 

We consider a random variable X belonging to the domain T>A and such that 
we are able to simulate X, T[X] and A[X]. 

Example 1. Wiener space. 

Let us consider as first example a stochastic differential equation defined on 
the Wiener space equipped with the Onstein-Uhlenbeck structure (cf [12], [4], [5]) 

(1) X t = x + f a(X s , s)dB s + f r(X s , s)ds 

Jo Jo 

By the functional calculus for the operators T and A (cf [7], [4]), if the coefficients 



are smooth, the triplet (X t , r[X t ], A[X t ]) is a diffusion solution to the equation 
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Denoting Y t the column vector (X t , T[X t ], A[X t ]) this equation writes 

(2) Y t = Y + f a(Y s , s)dB s + /* b(Y„ s)ds 

Jo Jo 

and solving it by the Euler scheme with mesh ^ on [0,T], i.e. 

[ns] 



(3) 



Y t n = Y + f a (Y£ s] , [ ^-)dB s + f b(Yr n su—)ds 
Jo ~~n~ n Jo — n 



yields a process Y t n = (X™, (r[X])™, (A[X])™)'. Now, it is straightforward to verify 
that the second and third components of this process are respectively equal to 
r[X f n ] and A[X™]. In other words, if for a process Z solution to a stochastic 
differential equation, we denote Z n the solution to the discretized s.d.e. by the 
Euler scheme of mesh ^ sur [0,T], we may write 



(4) 



T[X n ] 
A[X n ] 
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Thus, in order to compute the density of Xt, if we approximate it by the Euler 
scheme X^ and use the fact that the densities pr(xo,x) and p^(xo,x) of Xt and 
Xj< are close together and satisfy under regular hypotheses 

sup \p T (x ,x) -p T (x ,x)\ < — 

xo,x n 

(cf [1], [2], [10], [11, thm 4.1]) we are eventually in a situation where we have to 
estimate the density of X% in a framework where we are able to simulate X%, 
r[X$\ and A[X$\ thanks to the relations (4). 

Remark. Starting from the same equation (1), instead of putting an error on the 
Brownian motion, we can simply put an error on the initial value Xq. We obtain 
that (X t , T[X t ], A[X t ]) is still a diffusion, evidently different from the preceding 
one, but relations (4) are still valid, so that we know to simulate X%, T[X^\ and 

a \ x t\- , ^ ' 

In fact, the sequel will show that we have interest to consider both an error on 

the Brownian motion and an independent error on the initial value because this 

increases V\X\. 

Example 2. Poisson space. 

Let (H d ,i3(IR d ),/i,rf,7)b e an error structure on IR d whose generator is denoted 
(a, Va) and let iVbea Poisson point process on JR d with intensity measure ji. The 



space of definition of N, (Q, A, P), may be equipped with a so-called "white" error 
structure A, IP, E>, T) (cf [4]) with the following properties 

\/h E Va 

N{h) E ID and V[N(h)] = N(-f[h}) 
N(h) E VA and A[N(h)\ = N(a[h\). 

Simulating iV(£) amounts to draw a finite (poissonian) number of independent vari- 
ables with law /i, so that we are in a situation where N(h), T[N(h)], and A[N(h)] 
are simulatable, the same for a regular functional F(N(hx), N(h 2 ), . . . , N(h k )). 

Example 3. Monte Carlo space. 

Let X be a random variable simulatable on the Monte Carlo space by an 
infinite number of calls to the random function. Let us group together the calls 
with respect to which the variable X is regular and those with respect to which it 
is irregular or discontinuoous (use of the rejection method, etc.) so that X may 
be written on the space ([0, if, B([0, if), dx 1 ") x ([0, if, B([0, if), dx N ) 

X = F(U ,U 1 ,...,U m ,...;V ,V 1 ,...,V n ,...) 

where the C/j's are the coordinates of the first factor and the V}'s those of the 
second one, the function F being regular with respect to the U^s. 
Let us put on the U^s the error structure 

([0, If, B([0, If), dx™, ID, T) = ([0, 1], B([0, 1]), dx, d, 7 ) N 

where (d, 7) is the closure of the operator 7[w](x) = x 2 (l — x) 2 u' 2 (x) for u G 
^([0,1]). 

If F is C 2 with respect to each C/j sur [0, 1] and if the series 

00 1 

E(9*OT(! " U t ) 2 + m(l -Ui)(l- 2170) 
i=0 z 

converges in L 2 , we have, X e T>A and 

4*] = £&,(i*OT(i - u i) 2 + - Ui)(l - 2^)) 

so that X, T[X] et A[X] are simulatable. 
2. DIMINISHING THE BIAS. 

Let be (fi, ^4, P, D, T) an error structure, (£, ID) the associated Dirichlet form, 
and (A, VA) the associated generator. 

Let us explain the intuitive way. The symmetric Markov process associated 
to the error structure, in short time e, induces an error on any regular random 
variable defined on (Q,A,TP) whose bias is eA[X] and whose variance is eT[X]. 

Since the probability IP is invariant by the transition semi-group of the Markov 
process, the law of the random variable X is nearly the same as that of 



X + eA[X] + ^^T[X]G 



where G is an exogenous reduced Gaussian variable independent of A. It follows 
first that the random variable X + eA[X] which has the same expectation as X, 
possesses a smaller variance than that of X. This is shown by the following result. 

For X G (VA) d , we denote var[X] the co variance matrix of X, A[X] the col- 
umn vector with components (AfXj, . . . , A[X rf ]), T[X] the matrix r[Xj,X,-] and 

y / TjX]' the positive symmetric matrix square root of T[X]. 
Proposition 1. For X e (VA) d , 

d 

tiace(^[X + eA[X]\) = trace (y^i[X]) - 2e]T£[X i ] + e 2 ||A[X]|| 2 . 

i=i 

If A[X] is not zero, this quantity is minimum for e = J2i £[Xi]/\\A[X] \\ 2 and is 
equal to 

trace( m [X])-2^||. 
Proof. The result comes directly from the relation £[Xj\ = — < A[Xi], Xi >. o 

In order to calculate IEX, it is interesting to simulate X + eA[X] instead of X 
as soon as e e]0, 2 Ei £[Xi]/\\A[X] || 2 [. 

We apply now the same idea to the computation of the density of X that we de- 
note / when it exists. Let g(x — m, S) be the density of the normal law on TR d with 
mean m and covariance matrix E supposed to be invertible. Given X, A[X],'£pC] 

the conditional law of the random variable X + + v / £y / £[X]"G where G is 

an independent reduced Gaussian variable, has a density g(x — X — eA[X], e£.[X]) . 
The goal is to show, under suitable hypotheses assuring T\X] to be invertible, that 
JEg(x — X — eA[X], e£[X]) converges to f(x) faster than in the classical kernel 
method. 

Lemma 1. Let X be in (VA) d . We suppose that X possesses a conditional density 
r)(x, 7, a) given £[X] = 7 et A[X] — a such that x 1— > 77(2;, 7, a) be C 2 with bounded 
derivatives. Then Vrc e H d 

JE[-(A[X\yV xV (x,l[X],A[X\) + itrace (r[X] .Hess^) (x,L[X],A[X])] = 0. 

Proof. Let us give the argument in the case d — 1. Let be C 2 with compact 
support on JR. By [7], denoting A^ the generator in the L 1 sense, we have 

A^[p(X)] = v /(X)A[X\ + ^(X)T[X\. 

Hence if /x(c?7, da) is the law of the pair (T[X], 

IEL4 (1) [ip(X)\ = = J n{di, da) J (ip'(x)a + ip"(x)^)r](x, 7, a)dx. 

Integrating by parts gives 

J n{dry,da) J (p(x)(-arf x (x, 7, a) + \pr(' x i{x, 7, a))dx = 



hence ^ 

J fi(d'y,da)(-arj' x (x,r/,a) + -777^ (x, 7, a)) = 

as soon as E| -A[X]rj' x (X, F[X],A[X]) + \r[X]rf^(X, T[X], A[X})\ e L} oc (dx) what 
is satisfied under the assumptions of the statement. o 

First we study the bias : 

Proposition 2. Let X be as in the above lemma and let the conditional density 
r)(x, 7, a) be C 3 bounded with bounded derivatives. 
As e — > 0, the quantity 

j 2 (m\g{x -X- eA[X],eL[X])] - f(x)) 

possesses a finite limit equal to 

lJEU[X]y(Ress xV )(x,L[X},A[X])A[X] 

-E^ =1 ^[^]r[x,,x fc ]<^ fc (a;,r:[x],A[x])]. 

Demonstration. The argument begins with the relations 

JE[g(x - X - sA[X] , e£[X])] = / ^(cfry, da) J g{x - y - e, £7)77(2/, 7, a)dt 
= J //(c?7, da)JEr](x — ea — y/e^fyG, 7, a) 

where G is a reduced Gaussian variable with values in H rf and then consists of 
expanding i](x — ea — yfe^G, 7, a) by the Taylor formula and taking the expec- 
tation. The term of order zero gives f(x), the term in y/e vanishes since G is 
centered, the term in e is zero because of lemma 1, the term in e^fe vanishes be- 
cause G 3 is centered. The hypotheses give the upper bounds allowing to conclude. o 

About the variance we have : 

Proposition 3. Let X be as in proposition 1 and such that (det£[X]) _ ^ e L 1 , 

\rni e ^E d l 2 Mg 2 {x - X - eA[X},eZ[X\) = \im £ ^ e d/2 ^g(x - X - eA[X],eZ[X\) 

= V (x,g[X],A[X]) 

(47r) d / 2 y / det£[x] 
Demonstration. We have 

1Eg 2 (x- X -eA[X],eL[X]) = J pi(d-f,da) J g 2 (x - y - ea, £7)77(2/, 7, a)dy. 

Since ^ 

g2[ ~ Z ^ = (27ry/*(2ey/ 2 ^^ 9{z, 2 l) 
we obtain the result by dominated convergence and the continuity of r\. o 



3. COMPARISON OF RATES. 



The quantity JEg(x—X—eA[X],e^_[X]) is computed by the law of large numbers 
so that the approximation f(x) of f(x) is 

1 N 

/(*) =vE^- X n ~ <A[X]) n , e(L[X]) n ) (*) 

JV n=l 

where the indices n denote independent drawings. 

• If we are using the L 2 criterion 

ll/>)-/Wlli 2 =var/>) + (bias) 2 
we are led to choose e = \ and 

||/(x)-/(x)|| La = -i-0(l) 
to be compared with l 2 = 1 4 in the case of the classical kernel method (cf 

Nd+A J\f2d+8 

[13] [14]). We see that the new method divides the dimension by two. 

• The other criterion 

c(f(x), f(x)) = sup |IEV(/>) " ¥>(/(*))! 

where V is the set of polynomials of second degree (p(x) = ax 2 + bx + c with \a\ < 1 
and |6| < 1, what gives 

C (f(x)j(x)) = \nf 2 (x)] - + - mi 

may be better adapted to the case of error calculus for the reason that, when the 
errors are thought as germs, in short time, of Ito processes, biases have the same 
order of magnitude as variances (not as standard deviations). This criterion leads 
us to take e = — ^— what gives c(f(x),f(x)) = — h^Oil) to be compared with 

\ = 1 4 in the classical case, we see that for this criterion too the proposed 

JVd+2 J\r2ci+4 

method divides the dimension by two. 
4. DIRECT FORMULAE 

When X is real valued, explicit formulae may be proved that allow simulations 
at the speed of the law of large numbers, provided that in addition to X, T[X], 
and A [X], we are able to simulate the random variable T[X, fpq], what is possible 
under additional regularity assumptions. 

For instance in the example 2, we have easily, if X = N(h) 

1 N^jh^jh]]) 

[ T[x] J (AWD) 4 ■ 



Proposition 4. a) If X is in DA with T[X] E ID and T[X] > a.s. then X has 
a density f possessing an l.s.c. version f and 

(5) hx)=\im^v(sign(x-X)(T[X,—±— T ] + 2A[X] ^ 



eio 1 2 V 7V L 'e + T[X] J e + r[X]7' 

6j J/m addition £ ID, t/ien X /ias a density f which is absolutely continuous 
and 

(6) /(*) = (sign(a; - X)(T[X, J^] + ^S)) . 

Demonstration. Let us begin with the case a). Since X e "PA and T[X] e L 2 , 
for any C 2 function with bounded derivatives (cf [7] chap I), we have <p[X] G VA 
and 

AM*)] = <f/(X)A[X] + \v"{X)T[X] 

hence Vt > 

(7) * (x) = • 

Since E 2 ff|ffi = -JET[<p(X), ^pj] taking the expectation we obtain 



IX = 



Let us put K £ {x) = n^§q\ X = *] and H e (x) = E[(L[X, ^] + 
x]. Relation (5) writes 

(9) | ^"(:r)X £ (:r)IP x (d:r) = - | ^'(^ £ (x)P x ((i4 

The derivative in the distributions sense of the measure K £ (x)W x {dx) is the mea- 
sure H £ (x)W x (dx). It follows that the measure K £ (x)W x (dx) has a density and 
since K £ > Pj-a.s. (because T[X] > P-a.s.) the law Fx has a density /. 

(We prove here again the implication X e ID, T[X] > =>■ Px « dx which 
is true for any local Dirichlet form with square field operator cf [7]) 

Hence H £ (x)JP x (dx) = H £ (x)f(x)dx and (9) implies that the measure^ (x)JP x (dx) 
has an absolutely continuous density F £ {x) and Me > we have f(x) = ^ for 
almost every x. 

Taking ip(y) = yj A 2 + (y — x) 2 in (8), it comes 



V(A 2 + (X - x ) 2 )ie + T[X]J 



! " V =(r[*,— Ul + 2A[X] 



^x 2 + (x- x ) 2K 1 'e + nxy e + r[xy 



When A — > 0, by dominated convergence, for all x, the right-hand side converges 
to 

EM,-x)(r[x^H^)l 

where sign(y) = if y ^ and sign(0) = 0. 



The left-hand side is equal to 

/ w^rt K ' (v)Wx(dy) = I mlhrt FM dv 

since F £ is continuous, this converges when A — > to 2F £ (x), therefore we have 
the following equality between continuous functions 



Fc(x) = \e 



Now, as e [ 0, by its definition the function K £ increases to 1 Pi-a.s. since V[X\ is 
supposed to be strictly positive a.s. Hence K £ (x)f(x) increases to f(x) for almost 
every x. The equality F £ (x) = K £ (x)f(x) valid for almost every x implies that F £ 
is almost everywhere, hence everywhere, increasing and converges to / l.s.c. equal 
to / almost everywhere. 

In order to prove the point b) we proceed similarly and the hypotheses allow 
to replace (8) by the relation 

(10) E[/(I)] = -nv'(X)(T[X, -L] + ^E)]. 

Putting H(x) = E[r[X, + j^-\X = x] we see that the law of X, V x {dx), 
possesses a derivative in the sense of distributions H(x)JPx(dx) which is absolutely 
continuous, henceX has an absolutely continuous density /. Taking again ip(y) = 



X 2 + (y — x) 2 , we obtain 



by the same argument as above. o 

The density of X being obtained, we can extend the formulae (5) and (6) in 
order to compute conditional expectations. 

Proposition 5. Let be G G ID n L°° , 

a) under the assumptions of prop. 4 a ), we have dx-a.e. 

(11) f( x) nG\X = *] = lim i E (sign(, - X) (m -g^l + f^M)) . 
the right-hand side is l.s.c. if G > 0, 

b) under the assumptions of prop. 4 b) and with £ DPI L°° , we have dx-a.e. 

(12) f(xMG\X = x] = ilE (sign(x - X)(F[X, ^] + ^^)) ■ 
where the right-hand side is continuous. 



The proof is similar to that of proposition 4. 



o 



Let us remark eventually that letting tp'(X) go to 1 in formulae (8) and (10) 
and in the analoguous formulae of proposition 5 we see that VG e ID n L°° 



E(r[x, 



^ + r[x] 



G 



] + 



2GA[X] 

£ + r[x] 



) = o 



and also for e = under the hypotheses of prop. 4 b). Hence it is possible to 
introduce, as remarked in [10], an arbitrary control deterministic function c in 
order to optimize the variance. For instance formula (12) becomes 



Comment. In the kernel method (cf [13] [12] [8]), cancelling the first term in the 
asymptotic expansion of the bias is an old idea and has been explored by several 
authors either by the use of non-positive kernels (cf [9] [13]) either by a Romberg 
method what amounts to the preceding case. In the method we propose in sections 
2 and 3, the kernel is random and depends on the random variable itself. That 
shifts from an order of magnitude. Then the above idea may be applied again. 

The nearest work to the section 4 is the study by A. Kohatsu-Higa and R. 
Pettersson [10] which uses integration by parts on the Wiener space in the sense of 
Malliavin, also the paper of Bouchard, Ekeland and Touzi [3] . The difference in the 
points of view comes mainly from the fact that the integration by parts formulae 
are not the same, ours are simpler and do not involve Skorokhod integrals. 

Let us quote also that the results of sections 2 and 4 may be applied as well to 
the Poisson space or the Monte Carlo space with a possible choice of the Dirichlet 
form what gives a usefull flexibility in order to take in account the studied specific 
model. 
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